Theoretical Insights into the Ultrafast Deactivation Mechanism and Photostability of a Natural Sunscreen System: Mycosporine Glycine

In this work, different levels of quantum computational models such as MP2, ADC(2), CASSCF/CASPT2, and DFT/TD-DFT have been employed to investigate the photophysics and photostability of a mycosporine system, mycosporine glycine (MyG). First of all, a molecular mechanics approach based on the Monte Carlo conformational search has been employed to investigate the possible geometry structures of MyG. Then, comprehensive studies on the electronic excited states and deactivation mechanism have been conducted on the most stable conformer. The first optically bright electronic transition responsible for the UV absorption of MyG has been assigned as the S2 (1ππ*) owing to the large oscillator strength (0.450). The first excited electronic state (S1) has been assigned as an optically dark (1nπ*) state. From the nonadiabatic dynamics simulation model, we propose that the initial population in the S2 (1ππ*) state transfers to the S1 state in under 100 fs, through an S2/S1 conical intersection (CI). The barrierless S1 potential energy curves then drive the excited system to the S1/S0 CI. This latter CI provides a significant route for ultrafast deactivation of the system to the ground state via internal conversion.


INTRODUCTION
Ultraviolet radiation (UVR) is the most energetic region of the broad spectrum of wavelengths from solar radiation that reaches the earth. This radiation is further subdivided into UVA (400−315 nm), UVB (315−280 nm), and UVC (280− 100 nm). 1,2 Almost all the UVC and a large fraction of the UVB radiation are absorbed by the ozone layer in the stratosphere. This results in 5% UVB and 95% UVA accounting for the total UVR that reaches the earth's surface. These energy components have an impact on the earth's biosphere. 3 Although UVR can be beneficial, for example, the production of vitamin D being essential for prevention of skeletal disease, 4−6 there are lots of concerns regarding its hazardous effects on life. 7−10 The damaging effects of overexposure to UVR have been widely reported in previous studies and reviews, including cataract formation, skin aging, DNA mutation, and skin cancer. 11−15 Hence, the need for a balance between exposure to UVR and protection against overexposure to UVA and UVB is of crucial importance. 16 There are different types of photoprotective compounds present on the market. Nevertheless, so far, no compound has proved to be exceptionally versatile. A major concern regarding photoprotective systems (sunscreens) is that some of the active ingredients within sunscreens are inherently not photostable. A well-known example of this issue is the UV filter avobenzone. 17 A good candidate for an efficient UV-filter (the active ingredient within sunscreens) should have efficient nonradiative decay and long-term photostability to dissipate the excess absorbed energy safely and quickly to regenerate the original ground-state molecule ready for reabsorption of another photon. 12 Photoprotection in micro-and macroorganisms is achieved by a number of different molecules. 18 The particular family of molecules of interest here are mycosporines and mycosporinelike amino acids (MAAs). Mycosporines are derived from a cyclohexenone core, and MAAs are derived from a cyclohexenimine core. 19 Currently, there are over 70 mycosporines and MAAs isolated, and the various substituents on the ring, for example, the amino and imino functionalities, which are often amino acids or amino alcohols, are responsible for the wide range in peak absorption wavelength across both UVA and UVB regions of the electromagnetic spectrum. 20 These properties ideally lend themselves toward use in future cosmeceutical applications, largely due to the ability to tune the range (UVA and UVB) over which particular mycosporines and MAAs absorb. In contrast to the crucial photoprotection role that they have, rare reports in the literature are devoted to the fundamental root of the photophysical nature of these systems especially involving the cyclohexenone core. 21−23 Recent studies on MAAs clearly show that the excited electronic states are very short-lived, decaying rapidly to repopulate the ground electronic state. However understanding of the mechanisms of the relaxation dynamics in mycosporines and MAAs is still in its infancy. 24,25 Herein, we report on a systematic investigation of the photochemical mechanisms for radiationless excited-state deactivation of the most stable conformer of mycosporineglycine (abbreviated as MyG henceforth). This system has an absorption maximum at 310 nm in water solvent medium 26 and possesses a cyclohexenone ring linked with an amino acid, glycine ( Figure 1). However, to the best of our knowledge, there are no reports on the spectroscopy of this system in the gas phase. The motivation behind such studies is that it allows us to distinguish between intrinsic molecular properties and properties resulting from the biological environment. 23,27 The details of fundamental biochemical characteristics and interactions are often hidden by macroscopic solvent effects and interactions with other molecules. Thus, studying individual biomolecules in vacuum reveals the intrinsic properties of the most basic biomolecular processes. 27,28 In the present study, we explored the relevant excited-state potential-energy functions along a photochemical reaction path and optimized the relevant conical intersections (CIs) involved in the deactivation process. We have performed these explorations with the ADC(2) and CASSCF electronicstructure methods. In addition to the static exploration of the excited-state potential-energy profiles, we report the results of nonadiabatic trajectory-surface-hopping molecular dynamics simulations of MyG. Taking the implicit ethanol solvent model into account, the dynamics simulations were performed with the TD-DFT electronic-structure method, which is the method of choice for systems of this size.

COMPUTATIONAL DETAILS
The conformational landscapes of MyG have been determined based on the Monte Carlo multiple minimum (MCMM) methodology 29 using the MacroModel suite of program (version 8.5). 29 This initial exploration of the PES was done using classical molecular mechanics force fields (MMFFs). 30 The optimized geometry of selected conformers of the ground state and corresponding vibrational frequencies have been determined either at the density functional theory (DFT)/ B3LYP or RI-MP2 level of theory using the cc-pVDZ basis set. 31 The vertical and adiabatic excitation energies of the lowest excited singlet states and oscillator strengths for electronic transitions were computed based on the RI-ADC(2), TD-DFT/ωB97XD, 32 and multistate complete active space second-order perturbation (MS-CASPT2) 33 theory. The RI-MP2 and RI-ADC(2) calculations were carried out by the Turbomole program suit (V 6.3.1). 34,35 The DFT and TD-DFT calculations were employed using the Gaussian 16 program, 36 and CASSCF/CASPT2 37,38 computations were performed using OpenMolcas (V 18.09). 39,40 We have determined the potential energy curves representing the deactivation mechanism of our selected system based on linear interpolation in internal coordinates (LIIC). Our method of choice in this case is ADC(2) because of its reliability from our previous works along with works from other groups. 41−46 In addition, CI geometries were optimized using a state-averaged complete active space self-consistent field (SA-CASSCF). 47  Moreover, determination of the UV absorption spectrum and nonadiabatic dynamics simulation have been conducted based on the TD-DFT model using Newton-X (V. 2.0) 48 interfaced with Gaussian 16. The validity of the selected model for determining the excited-state deactivation processes has been confirmed in different reports. 16,49−52 To perform nonadiabatic dynamics simulations, the initial conditions were computed for geometry and nuclear momenta sampled from a Wigner distribution based on the S 0 normal modes. Vertical excitation energies were computed at the ωB97XD/6-31G* level, and these energies were obtained with δ = 0.05 eV Lorentzian line broadening for the two lowest singlet excited states.
In addition, to investigate the solvent's influence on absorption spectra, the polarizable continuum model (PCM/ ethanol model) 53 implemented in Gaussian 16 has been used.

Ground-State Optimized
Geometries. The conformational flexibility of MyG mainly arises from (1) the −glycine side chain, (2) the −methoxy group (−OCH 3 ), (3) the hydroxymethyl group (−CH 2 OH), and (4) flexible carbon sites of the ring. The orientation in space of the glycine side chain with respect to the six-membered ring could be subject to diversity in conformation in addition to location of its own atoms with respect to each other. Over the cyclohexenone ring, there are two double bonds, C 5 �C 6 and C 4 �O 9 ; thus, the C 4 , C 5 , and C 6 locates mainly in the same plane owing to the sp 2 hybridization nature of the C atoms (see Figure 1 for numbering). Other carbon atoms in the ring (C 1 , C 2 , and C 3 ) are sp 3 -hybridized, and there are two possibilities for outof-plane location of the C 1 atom and its subgroups; it could be pointing upward or pointing downward of the ring. Consequently, this offers another root of obtaining new conformations.
The conformational landscape of MyG has been determined based on the conformation search algorithm using the MCMM methodology. MMFFs were used for 1000 steps. All structures generated with energies within 20.0 kJ mol −1 of the minimum were examined. Based on the MCMM theoretical model, 101 conformers have been obtained. We have selected 30 of the most stable structures and determined the minimum geometry and consequently the relative stabilities based on the DFT method (at the B3LYP/cc-pVDZ level). Optimized geometries of the selected lowest-lying conformers (in energy) are presented in Figure 2. As shown, we have called the most stable structure MyG-A, which is 3.57 kJ mol −1 more stable than the next nearest conformer (at the B3LYP/cc-pVDZ level of theory, taking the ZPE correction into account). Comprehensive information regarding the 30 conformers is also presented in Table S1, Supporting Information file.
As shown in Figure 2, MyG-A exhibits a puckering feature around C 1 . The dihedral angle of C 1 −C 2 −C 3 −C 4 describing this deformation from the ring has been predicted at 49°. Other carbon sites of the six-membered ring stay roughly on the same plane; thus, the system looks like a half-chair. In addition, the dihedral angle of C 15 −O 10 −C 5 −C 6 describing the position of the −O−CH 3 side chain with respect to the cyclohexenone ring has been predicted at 115°which is nearperpendicular (downward). The glycine side chain remains in the same plane with respect to the six-membered cyclohexenone ring. Moreover, from inspection of Figure 2, it is apparent that the out-of-plane deformation of C 1 above the molecular plane is present in several structures (for instance, A, C, D etc.), while its out-of-plane deformation below the molecular plane is also significantly populated (such as B, E, I, L, and N). Furthermore, when the glycine side chain remains in the molecular plane, the corresponding conformers are significantly more stable than those that have other orientations (such as K, L, M, and N).
To note, we have optimized the 6 most stable conformers based on the MP2/cc-pVDZ theoretical model, with their relative stabilities confirming the DFT level results. The optimized structure of the most stable MyG-A conformer at the MP2 level has been determined to be fairly similar to that found with DFT, with slight differences in the angles and The MS-CASPT2 results have been obtained based on 8 electrons and 8 orbitals as the active space and 5 multistates. b Adiabatic S 1 ← S 0 transition energy determined at the optimized S 1 geometry of MyG-A. Further information regarding the MOs and transition characters could be found in Table S3, Supporting Information file.  Figure S1 of the Supporting Information. From our conformational study, we have selected the most stable conformer (MyG-A) and concentrated on its electronic structure, excited states, and deactivation mechanism which now follows.

Electronic Transition Energies.
Calculated vertical transition energies and oscillator strengths of the most stable conformer, MyG-A, to the four lowest lying singlet excited states (S 1 −S 4 ) have been determined and are presented in Table 1. We have employed different levels of theory: TD-DFT using the ωB97XD functional; RI-ADC(2); and MS-CASPT2. Regarding the TD-DFT calculations, the ωB97XD functional was used as it successfully described the energy of the electronic states and nonradiative deactivation process in similar systems. 24,49,51,54−56 Selected valence molecular orbitals (MOs) of MyG-A are presented in Figure 3. According to our ADC(2) results, the S 1 ← S 0 electronic transition has been assigned as an optically dark ( 1 nπ*) transition, originating from LUMO ← HOMO − 1 (85%) single-electron transitions. The HOMO and LUMO, respectively, stand for the highest occupied molecular orbital and the lowest unoccupied molecular orbital. The second singlet electronic transition (the S 2 ← S 0 ) has been assigned as an optically bright 1 ππ* transition, involving a LUMO ← HOMO transition (92%). As the molecular point group is C1 (and there is no symmetry plane), our assignments are only qualitative.
The MOs involving the S 3 ← S 0 electronic transition have been assigned mainly as LUMO + 1 ← HOMO and LUMO + 2 ← HOMO transitions, and for S 4 ← S 0 , the LUMO ← HOMO − 3 transition has been predicted to be the main contributor. More information regarding the electronic transitions and contributing MOs can be found in the Supporting Information.
The MS-CASPT2 results have been obtained, selecting the 8 electrons and 8 valance orbitals (4 occupied and 4 virtual) as the active space (see Table S2, Supporting Information). We have examined different active spaces to determine the vertical transition energies. As the S 2 −S 0 , 1 ππ* electronic transition dominates the photophysics of MyG, we have selected the active space yielding the most reliable result for this electronic transition. Although there is no experimental UV absorption spectrum of MyG in the gas phase, as mentioned before, the λ max of absorption around 300 nm has been reported for MyG in water solution. The active space of (8,8) gives better results than others, and no improvement has occurred considering the larger active spaces (see Table S4, Supporting Information file).
As shown in Table 1, the MS-CASPT2 results for S 1 and S 2 states are comparable with the ADC(2) and also TD-DFT/ ωB97XD results. There is a discrepancy relating to the S 3 ← S 0 electronic transition energy which is significantly underestimated at the TD-DFT level. However, due to the optically dark nature of this state (low oscillator strength), it should have no significant bearing on the photophysics of MyG. The energetic value of the S 2 ← S 0 transition ( 1 ππ*), responsible for the UV absorption of MyG-A, has been determined as 4.71, 4.69, and 4.68 eV, respectively, at the TD-ωB97XD, RI-ADC(2), and MS-CASPT2 levels of theory. Good agreement of TD-ωB97XD results with RI-ADC(2) and MS-CASPT2 (especially for the S 1 and S 2 transition energies) indicates that this functional is sufficiently reliable to describe electronic structures of the selected systems.
To perform a comparison between the five most stable structures of MyG-(A, B, C, D, and E), we have determined the vertical electronic transition energies in the gas phase and also with the implicit ethanol solvent for the four lowest lying singlet excited states at the TD-ωB97XD/cc-pVDZ level of theory. The results are presented in Table 2. Inspection of these results reflects at least two important points: (1) The electronic transition energies of the five most stable conformers of MyG are quite close to each other. This indicates that the photophysics of these systems should be similar. (2) The electronic transition energies in the implicit solvent (ethanol) are mildly different from the corresponding gas-phase results. Consistent with previous results, 57 it indicates that the solvent effect on the photophysical nature of these systems would be negligible.
To garner more insights into the exited-state nature of MyG-A, we have performed the RI-ADC(2) geometry optimizations for the excited S 1 and S 2 states. As mentioned supra, the S 1 is an optically dark state having 1 nπ* character and roughly zero oscillator strength. The S 2 ( 1 ππ*) state is responsible for the UV absorption of MyG (F = 0.450 in MyG-A). The S 1 excited state exhibits a local minimum, with a molecular geometry that is significantly different from that of the ground state (see Figure 4a,b). The glycine side chain exhibits an out-of-plane movement, and there is a slight distortion of the ring around Calculations were performed in the gas phase and in the PCM/ethanol implicit solvent, determined at the TD-ωB97XD/cc-pVDZ level of theory.
The Journal of Physical Chemistry A pubs.acs.org/JPCA Article the C 6 and C 3 regions. These changes in geometry stabilize the S 1 excited state, so that the adiabatic S 1 ← S 0 transition energy has been predicted to be significantly lower (2.86 eV) than its vertical analogue (4.03 eV). Furthermore, our calculations predict no local minimum on the S 2 potential energy surface, and the S 2 geometry optimization results in strong deformation of the six-membered ring from the C 6 carbon site and other side chains.

Photophysical Behaviors: CIs and Potential Energy Profiles.
To gain more insights into the excitedstate relaxation mechanism of MyG-A, we have searched for possible CIs. We have located two CIs for MyG-A based on the SA-CASSCF(6,6)/cc-pVDZ theoretical model (see Figure  5a,b). The active space in the SA-CASSCF calculations contains 6 electrons in 6 MOs (3 occupied and 3 virtual  MOs, Table S2, Supporting Information). The CI 1 represents the S 2 /S 1 potential energy crossing, and the CI 2 describes the minimum structure of the S 1 /S 0 curve crossing. As shown in Figure 5, both optimized CI structures exhibit puckering of the ring at C 6 , in addition to out-of-plane movement of the glycine moiety. This deformation is more pronounced in CI 2 . Moreover, there is a significant alteration in the C 8 −O 7 − C 5 −C 6 dihedral angle (from 115.6°in the ground state to −113.7°). On the other hand, in CI 2 , the methoxy (−OCH 3 ) group moves toward the glycine side chain (upward the ring). This alteration is accompanied by a slight out-of-plane distortion of C 5 downward the ring. The optimized Cartesian coordinates of the CIs are presented in the Supporting Information.
It is well known that locating CIs does not warrant the nonradiative deactivation process since the presence of a large barrier could prevent this deactivation pathway. Thus, we have determined potential energy (PE) profiles for the ground and two electronic excited states (S 1 and S 2 ) based on the LIIC connecting the Franck−Condon (FC) region of MyG-A to CI 2 . The results are presented in Figure 5. As shown, the S 1 and S 2 potential energy profiles decrease in energy along the LIIC coordinate. The S 2 potential energy curve first crosses with the S 1 potential energy curve (obtaining the CI 1 in the multidimensional picture).
This trend continues until the end of the reaction coordinate. The ground-state PE profile increases in energy along the LIIC coordinate, ultimately crossing the S 1 PE profile. As shown, the S 1 PE curve starting from the FC region passes through a shallow local minimum and then increases to cross with the S 2 PE curve. This local minimum in the LIIC potential energy profile is in accordance with the results of obtaining a minimum geometry for the S 1 state as discussed supra. The S 1 /S 0 curve crossing at the end of the reaction coordinate provides an essential root for ultrafast deactivation of the excited state based on internal conversion to the groundstate S 0 potential energy surface. Since the predicted CI (around 3.0 eV) is located significantly lower than the S 2 vertical transition energy [of 4.69 eV at the ADC(2) level], one could envisage that excited-state population in S 2 will approach the CI especially since the S 2 PE profile is purely repulsive.
One may be concerned that ADC(2), as a single-reference method, would not be reliable enough to describe the photophysical nature, especially the S 1 /S 0 CI of MyG-A. To address this concern, we have determined the ground-and excited-state potential energy profiles of MyG-A along the  In order to shed more light into the relaxation dynamics of MyG, we have employed the nonadiabatic dynamics simulation based on the TD-DFT theoretical model on the most stable conformer, MyG-A. In this regard, a nonadiabatic surface hopping dynamics simulation was performed from the second excited-state S 2 ( 1 ππ*) at the TD-ωB97xD/6-31G* level of theory, with ethanol as the implicit solvent.
The photoabsorption spectrum of MyG-A, shown in Figure  6, was simulated using the nuclear ensemble approach. 61 The velocity Verlet 62 algorithm with a time step of 0.5 fs was used for the integration of Newton's equations. Since the first bright state in MyG-A is the S 2 ( 1 ππ*) state, the nonadiabatic dynamics simulation has been started from S 2 , using the decoherence-corrected 63 fewest switches surface hopping (DC-FSSH) approach. 64 A total of 52 trajectories starting from the optically bright S 2 state were considered. The initial conditions for the dynamics were sampled within the 4.55 ± 0.10 eV spectral window highlighted in Figure 6. We considered that a trajectory went back to the ground-state S 0 when the S 1 −S 0 energetic gap approached a threshold of 0.25 eV.
Surface hopping dynamics was simulated for a maximum of 1000 fs with a nuclear time step of 0.50 fs. 100% of optically bright S 2 ( 1 ππ*)-state population is transferred to the S 1 ( 1 nπ*) state within the initial ∼120 fs of the dynamics. Therefore, the deactivation mechanism starts at the S 2 ( 1 ππ*) state, with excited-state population traversing the S 2 /S 1 CI in around 120 fs, leading the excited system to drive over the S 1 PE curve ending by CI 2 (see Figure 6b). This later CI provides an essential root for ultrafast internal conversion of the excited system to the ground state. Consequently, the net effect of this process would be converting the hazardous UV light (of 4.69 eV) to harmless vibrational motion (or heat). The dynamics results are in good agreement with our previously discussed ab initio results.
From the results of the dynamics, we predict that 82% of the trajectories relax based on the ring puckering from the C 6 region (corresponding to CI 2 presented in Figure 5b) as well as the C 3 region. Moreover, 10% trajectories have been assigned to relax with ring puckering from C 1 and C 4 carbon sites.
To conclude this section, we present the time evolution of the nonadiabatic population of the ground and the first two excited states from the nonadiabatic dynamics simulation of MyG-A ( Figures S3 and S4, Supporting Information). Using a sigmoid curve fitting procedure 56 as an approximate guide, we determine that the S 2 excited-state lifetime; τ Sd 2 is predicted to be ∼67 fs, while the S 1 excited-state lifetime, τ Sd 1 , is predicted to be ∼517 fs. Further details regarding this fitting is presented in Supporting Information, Figure S2. From this, we conclude that photoexcitation of MyG to S 2 can relax to the S 0 ground state in approximately τ Sd 2 + τ Sd 1 = 584 fs. This ultra-short lifetime of the excited state is in accordance with the high photostability of MyG-A.

SUMMARY AND CONCLUSIONS
Molecular mechanics configurational search has been employed to determine the possible structures and conformers of mycosporine glycine. We have employed stepwise computational models to determine the lowest lying conformers and consequently the most stable structure. Following this, state-ofthe-art ab initio computational models as well as nonadiabatic dynamics simulations have been conducted to shed light on the photophysics of MyG. We show that the S 2 ( 1 ππ*) state is responsible for the UV absorption. Following photoabsorption in the UV region [4.69 eV at the ADC(2) theoretical level], excited population proceeds along a ring twisting reaction coordinate, traversing through the S 2 /S 1 CI. This takes place on the order of 100 fs. We have also shown that the potential energy profile of the S 1 state after the S 2 /S 1 CI significantly decreases in energy and thus drives the excited system to the S 1 /S 0 CI. This CI governs the ultrafast deactivation of the excited system to the ground state by internal conversion. These results provide significant molecular level insights into the origins of photostability of MyG and consequently its ability to photoprotect microorganisms against the potentially detrimental effects of UV radiation exposure. These results also provide benchmark data that could be used to develop nextgeneration UV filters for commercial applications.